####################################
###   REPLICATION MATERIALS FOR  ###
### EXECUTIVE COALITION BUILDING ###
####################################

#Load data#
setwd("/Users/napolio/Library/CloudStorage/Dropbox/Interagency Rulemaking/Analysis")
load("panel.RData")
load("panel_pres.RData")
load("panel_pres_slim.RData")
load("centrality.RData")
load("coalitions.RData")
load("oira_requests.RData")
load("reg_review.RData")


#Load packages#
library(ggplot2)
library(igraph)
library(stargazer)
library(plyr)
library(dplyr)
library(sandwich)
library(extrafont)
library(sampleSelection)

devtools::source_url("https://raw.githubusercontent.com/jkarreth/JKmisc/master/ggintfun.R")
ggintfun <- function (obj, varnames, varlabs, 
                      title = FALSE, 
                      rug = TRUE, 
                      twoways = TRUE, 
                      rugsize = 0.1, 
                      jitter_factor = 0, varcov){
  
  require(ggplot2); require(gridExtra)
  
  if(class(obj) == "lm"){
    
    if (!("model" %in% names(obj))) {
      obj <- update(obj, model = T)
    }
    
    v1 <- varnames[1]
    v2 <- varnames[2]
    vlab1 <- varlabs[1]
    vlab2 <- varlabs[2]
    
    if(length(unique(model.matrix(obj)[, paste(v1)])) == 2){
      s1 <- c(min(model.matrix(obj)[, paste(v1)]), max(model.matrix(obj)[, paste(v1)]))
    }
    
    if(length(unique(model.matrix(obj)[, paste(v1)])) > 2){
      s1 <- seq(from = min(model.matrix(obj)[, paste(v1)]), to = max(model.matrix(obj)[, paste(v1)]), length.out = 25)
    }
    
    if(length(unique(model.matrix(obj)[, paste(v2)])) == 2){
      s2 <- c(min(model.matrix(obj)[, paste(v2)]), max(model.matrix(obj)[, paste(v2)]))
    }
    
    if(length(unique(model.matrix(obj)[, paste(v2)])) > 2){
      s2 <- seq(from = min(model.matrix(obj)[, v2]), to = max(model.matrix(obj)[, v2]), length.out = 25)
    }
    
    
    
    b1.pos <- grep(pattern = v1, x = names(coef(obj)), fixed = TRUE)[1]
    b3.pos <- grep(pattern = v1, x = names(coef(obj)), fixed = TRUE)[2]
    b1 <- as.numeric(obj$coef[b1.pos])
    b3 <- as.numeric(obj$coef[b3.pos])
    var1 <- varcov[b1.pos, b1.pos]
    var3 <- varcov[b3.pos, b3.pos]
    cov13 <- varcov[b1.pos, b3.pos]
    
    eff1 <- b1 + b3 * s2
    var.eff1 <- var1 + s2^2 * var3 + 2 * s2 * cov13
    se.eff1 <- sqrt(var.eff1)
    low1 <- eff1 - 1.96 * se.eff1
    up1 <- eff1 + 1.96 * se.eff1
    
    b2.pos <- grep(pattern = v2, x = names(coef(obj)), fixed = TRUE)[1]
    b3.pos <- grep(pattern = v2, x = names(coef(obj)), fixed = TRUE)[2]
    b2 <- as.numeric(obj$coef[b2.pos])
    b3 <- as.numeric(obj$coef[b3.pos])
    var2 <- varcov[b2.pos, b2.pos]
    var3 <- varcov[b3.pos, b3.pos]
    cov23 <- varcov[b2.pos, b3.pos]
    
    eff2 <- b2 + b3 * s1
    var.eff2 <- var2 + s1^2 * var3 + 2 * s1 * cov23
    se.eff2 <- sqrt(var.eff2)
    low2 <- eff2 - 1.96 * se.eff2
    up2 <- eff2 + 1.96 * se.eff2
    
    rug2.dat <- data.frame(obj$model[v2])
    rug1.dat <- data.frame(obj$model[v1])
    
  }
  
  if(class(obj) == "lmerMod"){
    
    # Model matrix always included as obj@frame
    
    v1 <- varnames[1]
    v2 <- varnames[2]
    vlab1 <- varlabs[1]
    vlab2 <- varlabs[2]
    
    b1.pos <- grep(v1, names(fixef(obj)))[1]
    b3.pos <- grep(v1, names(fixef(obj)))[2]
    b1 <- as.numeric(fixef(obj)[b1.pos])
    b3 <- as.numeric(fixef(obj)[b3.pos])
    var1 <- varcov[b1.pos, b1.pos]
    var3 <- varcov[b3.pos, b3.pos]
    cov13 <- varcov[b1.pos, b3.pos]
    
    if(length(unique(model.matrix(obj)[, paste(v1)])) == 2){
      s1 <- c(min(model.matrix(obj)[, paste(v1)]), max(model.matrix(obj)[, paste(v1)]))
    }
    
    if(length(unique(model.matrix(obj)[, paste(v1)])) > 2){
      s1 <- seq(from = min(model.matrix(obj)[, paste(v1)]), to = max(model.matrix(obj)[, paste(v1)]), length.out = 25)
    }
    
    if(length(unique(model.matrix(obj)[, paste(v2)])) == 2){
      s2 <- c(min(model.matrix(obj)[, paste(v2)]), max(model.matrix(obj)[, paste(v2)]))
    }
    
    if(length(unique(model.matrix(obj)[, paste(v2)])) > 2){
      s2 <- seq(from = min(model.matrix(obj)[, v2]), to = max(model.matrix(obj)[, v2]), length.out = 25)
    }
    
    eff1 <- b1 + b3 * s2
    var.eff1 <- var1 + s2^2 * var3 + 2 * s2 * cov13
    se.eff1 <- sqrt(var.eff1)
    low1 <- eff1 - 1.96 * se.eff1
    up1 <- eff1 + 1.96 * se.eff1
    
    b2.pos <- grep(v2, names(fixef(obj)))[1]
    b3.pos <- grep(v2, names(fixef(obj)))[2]
    b2 <- as.numeric(fixef(obj)[b2.pos])
    b3 <- as.numeric(fixef(obj)[b3.pos])
    var2 <- varcov[b2.pos, b2.pos]
    var3 <- varcov[b3.pos, b3.pos]
    cov23 <- varcov[b2.pos, b3.pos]
    
    eff2 <- b2 + b3 * s1
    var.eff2 <- var2 + s1^2 * var3 + 2 * s1 * cov23
    se.eff2 <- sqrt(var.eff2)
    low2 <- eff2 - 1.96 * se.eff2
    up2 <- eff2 + 1.96 * se.eff2
    
    rug2.dat <- data.frame(obj@frame[v2])
    rug1.dat <- data.frame(obj@frame[v1])
    
  }
  
  plot1.dat <- data.frame(s2, eff1, se.eff1, low1, up1)
  
  plot2.dat <- data.frame(s1, eff2, se.eff2, low2, up2)  
  
  if (twoways == FALSE) {
    
    if(length(unique(model.matrix(obj)[, paste(v2)])) == 2){
      
      p1 <- ggplot(data = plot1.dat, aes(x = factor(s2), y = eff1), environment = environment()) +     # env. important for rug
        geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 0.5) + 
        geom_segment(aes(x = factor(s2), xend = factor(s2), y = low1, yend = up1), color = "black") +
        geom_point() + 
        xlab(paste(vlab2)) + ylab(" ")
    }
    
    if(length(unique(model.matrix(obj)[, paste(v2)])) > 2){
      
      p1 <- ggplot(data = plot1.dat, aes(x = s2, y = eff1), environment = environment()) +     # env. important for rug
        geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 0.5) + 
        geom_ribbon(aes(x = s2, ymin = low1, ymax = up1), alpha = 0.25, color = NA) +
        geom_line() + 
        xlab(paste(vlab2)) + ylab(" ")
    }
    
    if(title == TRUE){
      p1 <- p1 + ggtitle(paste("Marginal Effect of ", vlab1, sep = ""))
    }  
    
    if(rug == TRUE & length(unique(model.matrix(obj)[, paste(v2)])) > 2){
      p1 <- p1 + geom_rug(data = rug2.dat, aes(x = jitter(rug2.dat[, 1], factor = jitter_factor), y = 0), sides = "b", size = rugsize)
    }     
    
    p1 <- p1 + theme_bw()
    return(p1)
  }
  
  if (twoways == TRUE) {
    if(length(unique(model.matrix(obj)[, paste(v2)])) == 2){
      p1 <- ggplot(data = plot1.dat, aes(x = factor(s2), y = eff1), environment = environment()) +     # env. important for rug
        geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 0.5) + 
        geom_segment(aes(x = factor(s2), xend = factor(s2), y = low1, yend = up1), color = "black") +
        geom_point() + 
        xlab(paste(vlab2)) + ylab(" ")
    }
    
    if(length(unique(model.matrix(obj)[, paste(v2)])) > 2){
      p1 <- ggplot(data = plot1.dat, aes(x = s2, y = eff1), environment = environment()) +     # env. important for rug
        geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 0.5) + 
        geom_ribbon(aes(x = s2, ymin = low1, ymax = up1), alpha = 0.25, color = NA) +
        geom_line() + 
        xlab(paste(vlab2)) + ylab(" ")
    }   
    
    p1 <- p1 + ggtitle(paste("Marginal Effect of ", vlab1, sep = "")) + theme_bw()
    
    if(rug == TRUE & length(unique(model.matrix(obj)[, paste(v2)])) > 2){
      p1 <- p1 + geom_rug(data = rug2.dat, aes(x = jitter(rug2.dat[, 1], factor = jitter_factor), y = 0), sides = "b", size = rugsize)
    }     
    
    if(length(unique(model.matrix(obj)[, paste(v1)])) == 2){
      p2 <- ggplot(data = plot2.dat, aes(x = factor(s1), y = eff2), environment = environment()) +     # env. important for rug
        geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 0.5) + 
        geom_segment(aes(x = factor(s1), xend = factor(s1), y = low2, yend = up2), color = "black") +
        geom_point() + 
        xlab(paste(vlab1)) + ylab(" ")
    }
    
    if(length(unique(model.matrix(obj)[, paste(v1)])) > 2){
      p2 <- ggplot(data = plot2.dat, aes(x = s1, y = eff2), environment = environment()) +     # env. important for rug
        geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 0.5) + 
        geom_ribbon(aes(x = s1, ymin = low2, ymax = up2), alpha = 0.25, color = NA) +
        geom_line() + 
        xlab(paste(vlab1)) + ylab(" ")
      
    }  
    
    p2 <- p2 + ggtitle(paste("Marginal Effect of ", vlab2, sep = "")) + theme_bw()
    
    if(rug == TRUE & length(unique(model.matrix(obj)[, paste(v1)])) > 2){
      p2 <- p2 + geom_rug(data = rug1.dat, aes(x = jitter(rug1.dat[, 1], factor = jitter_factor), y = 0), sides = "b", size = rugsize)
    }     
    
    grid.arrange(p1, p2, ncol = 2)
  }
}

#Table 1#
coalitions[c(3,1,2,4,5),c(4,8,1,5,7,6)]

#Figure 1#
panel_pres_slim$dyad_type <- NA
panel_pres_slim$dyad_type[grepl("department", as.character(panel_pres_slim$agency_1), ignore.case = T) &
                            grepl("department", as.character(panel_pres_slim$agency_2), ignore.case = T)] <- "Cabinet-Cabinet"
panel_pres_slim$dyad_type[!grepl("department", as.character(panel_pres_slim$agency_1), ignore.case = T) &
                            !grepl("department", as.character(panel_pres_slim$agency_2), ignore.case = T)] <- "Independent-Independent"
panel_pres_slim$dyad_type[(grepl("department", as.character(panel_pres_slim$agency_1), ignore.case = T) &
                             !grepl("department", as.character(panel_pres_slim$agency_2), ignore.case = T)) |
                            (grepl("department", as.character(panel_pres_slim$agency_2), ignore.case = T) &
                               !grepl("department", as.character(panel_pres_slim$agency_1), ignore.case = T))] <- "Cabinet-Independent"
prop.table(table(panel_pres_slim$dyad_type))

by_type <- data.frame(rule = tapply(panel_pres_slim$rule, panel_pres_slim$dyad_type, mean))
by_type$dyad_type <- rownames(by_type)

ggplot(by_type) + geom_bar(aes(x=dyad_type, y=rule), stat="summary", fun.y=mean) + theme_bw() +
  xlab(" ") + ylab(" ") + ggtitle("Probability of Coalition by Agency Type Combination") + theme(panel.grid = element_blank()) +
  ylim(c(0,.6))

#Table 2#
panel_net <- subset(panel, dyad %in% panel_pres_slim$dyad)
net <- data.frame(weight = tapply(panel_net$weight, as.character(panel_net$dyad), sum))

net$dyad <- rownames(net)
net <- join(net, panel_net, by="dyad", type="left", match="first")
net$agency_1 <- as.character(net$agency_1)
net$agency_2 <- as.character(net$agency_2)
net2 <- net[,c(4,5,2)]
net2 <- subset(net2, weight>0)
g <- graph_from_data_frame(net2, directed=F)

g <- graph_from_data_frame(net2, directed=F)
head(sort(degree(g), decreasing = T), 5)
head(sort(degree(g), decreasing = F), 5)
head(sort(betweenness(g), decreasing = T), 5)
head(sort(betweenness(g), decreasing = F), 5)

#Figure 2#
ggplot(panel_pres_slim, aes(x=overlapping_laws, y=rule)) + 
  geom_point(pch=1) + 
  geom_smooth(method="glm", method.args=list(family="binomial"), color="black") +
  theme_minimal() + facet_wrap(~pres_term2, nrow=1) +
  xlab("Count of Laws Delegating to Both Agencies in Dyad") +
  ylab("Probability of Coalition Formation") +
  scale_y_continuous(breaks = c(0,.25,.5,.75,1),
                     labels = paste(c(0,25,50,75,100), "%"))

#Table 3#
panel_pres_slim$ag_align <- ifelse(panel_pres_slim$ideo_prox>mean(panel_pres_slim$ideo_prox), 1, 0)
panel_pres_slim$pres_align <- ifelse(panel_pres_slim$avg_pres_dist>mean(panel_pres_slim$avg_pres_dist), 0 , 1)
crosstab <- aggregate(rule ~ ag_align + pres_align, mean, data = panel_pres_slim)
crosstab$sd <- aggregate(rule ~ ag_align + pres_align, sd, data = panel_pres_slim)$rule
crosstab$n <- aggregate(rule ~ ag_align + pres_align, length, data = panel_pres_slim)$rule
crosstab$se <- crosstab$sd/sqrt(crosstab$n)

crosstab

#Table 4#
#Regress joint rule on main variables#
mod1 <- lm(rule ~ ideo_prox + avg_pres_dist + dyad + pres_term, data=panel_pres_slim)
mod2 <- lm(rule ~ ideo_prox*avg_pres_dist + dyad + pres_term, data=panel_pres_slim)
mod3 <- lm(rule ~ ideo_prox + avg_pres_dist + overlapping_laws + total_polit +  log(total_mentions+1) + log(total_rules_sum) + emp_diff + polit_diff + avg_house_dist + avg_senate_dist + avg_court_dist + dyad + pres_term, data=panel_pres_slim)
mod4 <- lm(rule ~ ideo_prox*avg_pres_dist + overlapping_laws  + total_polit +  log(total_mentions+1) + log(total_rules_sum) + emp_diff + polit_diff + avg_house_dist + avg_senate_dist + avg_court_dist + dyad + pres_term, data=panel_pres_slim)


#Caluculte cluster-robust standard errors by dyad#
vcov1 <- vcovCL(mod1, cluster=panel_pres_slim$dyad)
vcov2 <- vcovCL(mod2, cluster=panel_pres_slim$dyad)
vcov3 <- vcovCL(mod3, cluster=panel_pres_slim$dyad)
vcov4 <- vcovCL(mod4, cluster=panel_pres_slim$dyad)

ses <- list(vcov1, vcov2, vcov3, vcov4)
ses <- lapply(ses, diag)
ses <- lapply(ses, sqrt)

stargazer(mod1, mod2, mod3, mod4, se=ses, 
          keep=c("ideo_prox", "avg_pres_dist", "emp_diff", "polit_diff", "avg_house_dist", "avg_senate_dist",  "total_rules_sum",
                 "overlapping_laws", "total_mentions", "avg_court_dist", "total_polit"),
          dep.var.labels.include = T, keep.stat=c("adj.rsq", "n"), digits=3,
          star.cutoffs = c(.05, .01, .001), type = "text")

#Figure 3#
ggintfun(mod4, varnames=c("ideo_prox", "avg_pres_dist"), varlabs=c("Agency Alignment on Probability of Coalition Formation", "Presidential Misalignment \n (Centered on Zero)"),
         twoways=F, varcov = vcov4, title = T) +
  theme(panel.grid.minor = element_blank()) +
  scale_y_continuous(breaks=seq(-1,1,.25), name="Change in Probability of Coalition Formation\nGiven Unit Increase in Agency Alignment",
                     sec.axis = sec_axis(~. *.21, name = "Change in Probability of Coalition Formation\nGiven Standard Within-Dyad Increase in Agency Alignment",
                                         breaks=c(-.15,-.1,-.05,0,.05,.1,.15,.2))) +
  scale_x_continuous(breaks=c(-.5, -.25, 0,.25,.5)) + 
  geom_segment(aes(x=.308, xend=.308, y=-1, yend=.36),
               arrow=arrow(ends="last", angle=30, length=unit(.1, "in"))) +
  geom_label(x=.45, y=-.4, label="One Standard Deviation\nabove the Mean of\nPresidential Misalignment", family="Palatino Linotype") +
  geom_hline(yintercept=.4) + geom_label(family="Palatino Linotype", x=-.25, y=.625, label="Marginal Effect of Agency Alignment\nOne Standard Deviation above the Mean\nof Presidential Misalignment") +
  geom_point(aes(x=.308,.4), size=3, fill="white", pch=21)


#Table 5#

mod1 <- lm(change_request ~ coalition + as.factor(year), review_dat) 

mod2 <- lm(change_request ~ coalition + as.factor(year) + 
             ECONOMICALLY_SIGNIFICANT + MAJOR +
             FEDERALISM_IMPLICATIONS,
           review_dat)

se <- list(sqrt(diag(vcovCL(mod1, cluster = review_dat$year))),
           sqrt(diag(vcovCL(mod2, cluster = review_dat$year))))

stargazer(mod1, mod2, type = "text",
          keep = c("coalition", "ECONOMICALLY_SIGNIFICANT",
                   "MAJOR", "FEDERALISM_IMPLICATIONS"),
          se = se, keep.stat = c("adj.rsq", "n"),
          star.cutoffs = c(.05, .01, .001),
          dep.var.labels= "OIRA Requests Change",
          covariate.labels = c("Produced by Coalition",
                               "Economically Significant",
                               "Major", "Federalism Implications"))
#Table 6#
oira_request$cj_dist <- (oira_request$cj_dist-mean(oira_request$cj_dist))/sd(oira_request$cj_dist)

mod <- lm(coalition_change_request ~ cj_dist + politicization + 
            solo_change_request, oira_request)

mod2 <- lm(coalition_change_request ~ cj_dist + I(cj_dist^2) + politicization +
             solo_change_request, oira_request)


mod_sec <- heckit(coalition ~ cj_dist + I(cj_dist^2) + politicization + pres + name,
                  coalition_change_request ~cj_dist  + politicization + solo_change_request,
                  oira_request, method = "2step")


mod_sec2 <- heckit(coalition ~ cj_dist + politicization + pres + name,
                   coalition_change_request  ~cj_dist + I(cj_dist^2) + politicization + solo_change_request,
                   oira_request, method = "2step")

stargazer(mod, mod2, mod_sec, mod_sec2, type = "text",
          keep.stat = c("n", "adj.rsq"),
          keep = c("cj_dist", "politicization", "solo_change_request"),
          se = list(sqrt(diag(vcovCL(mod, oira_request$name))),
                    sqrt(diag(vcovCL(mod2, oira_request$name))),
                    sqrt(diag(vcov(mod_sec, part = "outcome"))),
                    sqrt(diag(vcov(mod_sec2, part = "outcome")))),
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels = "Pr(Change Request): Coalition Policy",
          covariate.labels = c("Presidential Misalignment", "Presidential Misalignment Sq.",
                               "Politicization", "Pr(Change Request): Solo Policy"))
#Figure 4#

predictions <- data.frame(cj_dist = rep(seq(-1.54,1.95,.01), 500))
predictions$fitted <- NA
for(i in 1:nrow(predictions)){
  b1 <- rnorm(1, mod_sec2$lm$coefficients[2], sqrt(diag(vcov(mod_sec2, part = "outcome")))[2])
  b2 <- rnorm(1, mod_sec2$lm$coefficients[3], sqrt(diag(vcov(mod_sec2, part = "outcome")))[3])
  b0 <- rnorm(1, mod_sec2$lm$coefficients[1], sqrt(diag(vcov(mod_sec2, part = "outcome")))[1])
  predictions$fitted[i] <- mod_sec2$lm$coefficients[1]*b0 + predictions$cj_dist[i]*b1 + 
    (predictions$cj_dist[i]^2)*b2 + mod_sec2$lm$coefficients[4]*mean(oira_request$politicization) +
    mod_sec2$coefficients[5]*mean(oira_request$solo_change_request, na.rm = T)
  if(i %% 1000 == 0){message(paste(i, nrow(predictions), sep = "/"))}
}

predictions_agg <- aggregate(fitted ~ cj_dist, mean, data = predictions)
predictions_agg2 <- aggregate(fitted ~ cj_dist, quantile, c(.05, .95), data = predictions)
predictions_agg2 <- predictions_agg2$fitted
predictions_agg <- cbind(predictions_agg, predictions_agg2)
names(predictions_agg) <- c("cj_dist", "y", "ymin", "ymax")

ggplot(predictions_agg, aes(x = cj_dist, y = y, ymin = ymin, ymax = ymax)) +
  geom_smooth(color="black") + geom_ribbon(alpha = .25) + theme_minimal() + 
  xlab("Presidential Misalignment\n(Centered on Zero)") +
  ylab("Predicted Probability of OIRA Review of Coalition Policy") +
  scale_x_continuous(breaks = seq(-1.5,2,.5)) + scale_y_continuous(breaks = seq(0,1.5,.25)) +
  geom_hline(yintercept = .5, lty = 2, col = "red", size = 1) +
  geom_rug(inherit.aes = FALSE, data = oira_request, aes(x = cj_dist), size = 2)


#APPENDIX#

#Table B1#
d_mod <- lm(degree ~ politicization + pres_prox + ln_emp + factor(year), centrality)
b_mod <- lm(betweenness ~ politicization + pres_prox + ln_emp + factor(year), centrality)
s_mod <- lm(I(log(strength+1)) ~ politicization + pres_prox + ln_emp + factor(year), centrality)

d_modFE <- lm(degree ~ politicization + pres_prox + ln_emp + factor(agency) + factor(year), centrality)
b_modFE <- lm(betweenness ~ politicization + pres_prox + ln_emp + factor(agency) + factor(year), centrality)
s_modFE <- lm(I(log(strength+1)) ~ politicization + pres_prox + ln_emp + factor(agency) + factor(year), centrality)

#Calculate cluster-robust standard errors by agency#

d_vcov <- vcovCL(d_mod, cluster=centrality$agency)
b_vcov <- vcovCL(b_mod, cluster=centrality$agency)
s_vcov <- vcovCL(s_mod, cluster=centrality$agency)

d_vcovFE <- vcovCL(d_modFE, cluster=centrality$agency)
b_vcovFE <- vcovCL(b_modFE, cluster=centrality$agency)
s_vcovFE <- vcovCL(s_modFE, cluster=centrality$agency)


d_se <- sqrt(diag(d_vcov))
b_se <- sqrt(diag(b_vcov))
s_se <- sqrt(diag(s_vcov))

d_seFE <- sqrt(diag(d_vcovFE))
b_seFE <- sqrt(diag(b_vcovFE))
s_seFE <- sqrt(diag(s_vcovFE))


stargazer(d_mod,d_modFE, b_mod, b_modFE, s_mod, s_modFE,
          se=list(d_se, d_seFE, b_se, b_seFE, s_se,s_seFE),
          keep=c("politicization", "pres_prox", "ln_emp"), nospace=T,
          digits = 3, digits.extra = 0, star.cutoffs = c(.05, .01, .001),
          keep.stat = c("n", "adj.rsq"),
          covariate.labels = c("Politicization", "Ideological Proximity to President",
                               "Logged Employees"), type = "text")

#Table C1#

max <- tapply(panel_pres_slim$rule, panel_pres_slim$dyad, max)
min <- tapply(panel_pres_slim$rule, panel_pres_slim$dyad, min)
range <- data.frame(max-min)
names(range) <- "range"
range$dyad <- rownames(range)
range <- subset(range, range!=0)

mod1l <- glm(rule ~ ideo_prox + avg_pres_dist + dyad + pres_term, data=panel_pres_slim[panel_pres_slim$dyad %in% range$dyad,], family="binomial")
mod2l <- glm(rule ~ ideo_prox*avg_pres_dist + ideo_prox + dyad + pres_term, data=panel_pres_slim[panel_pres_slim$dyad %in% range$dyad,], family="binomial")
mod3l <- glm(rule ~ ideo_prox + avg_pres_dist + avg_court_dist + overlapping_laws + log(total_mentions+1) + log(total_rules_sum) + emp_diff + polit_diff + total_polit + avg_house_dist + dyad + pres_term, data=panel_pres_slim[panel_pres_slim$dyad %in% range$dyad,], family="binomial")
mod4l <- glm(rule ~ ideo_prox*avg_pres_dist  + avg_court_dist + overlapping_laws + log(total_mentions+1) + log(total_rules_sum) + emp_diff + polit_diff + total_polit + avg_house_dist + dyad + pres_term, data=panel_pres_slim[panel_pres_slim$dyad %in% range$dyad,], family="binomial")


#Calculate cluster-robust standard errors by dyad#
vcov1l <- vcovCL(mod1l, cluster=panel_pres_slim[panel_pres_slim$dyad %in% range$dyad,]$dyad)
vcov2l <- vcovCL(mod2l, cluster=panel_pres_slim[panel_pres_slim$dyad %in% range$dyad,]$dyad)
vcov3l <- vcovCL(mod3l, cluster=panel_pres_slim[panel_pres_slim$dyad %in% range$dyad,]$dyad)
vcov4l <- vcovCL(mod4l, cluster=panel_pres_slim[panel_pres_slim$dyad %in% range$dyad,]$dyad)

ses <- list(vcov1l, vcov2l, vcov3l, vcov4l)
ses <- lapply(ses, diag)
ses <- lapply(ses, sqrt)

stargazer(mod1l, mod2l, mod3l, mod4l, se=ses,
          keep=c("ideo_prox", "avg_pres_dist", "emp_diff", "polit_diff", "avg_house_dist", "total_rules_sum",
                 "overlapping_laws", "total_mentions", "total_polit", "avg_court_dist"),
          dep.var.labels.include = T, keep.stat=c("ll", "n"), digits=3,
          star.cutoffs = c(.05, .01, .001), type = "text")


#Table C2#
mod1a <- lm(rule ~ ideo_prox + avg_pres_dist + dyad + pres_term, data=panel_pres)
mod2a <- lm(rule  ~ ideo_prox*avg_pres_dist + dyad + pres_term, data=panel_pres)
mod3a <- lm(rule  ~ ideo_prox + avg_pres_dist + avg_court_dist + overlapping_laws + log(total_mentions+1) + log(total_rules_sum+1) + emp_diff + avg_house_dist + dyad + pres_term, data=panel_pres)
mod4a <- lm(rule  ~ ideo_prox*avg_pres_dist + avg_court_dist + overlapping_laws + log(total_mentions+1) +log(total_rules_sum+1) + emp_diff + avg_house_dist + dyad + pres_term, data=panel_pres)


#Caluculte cluster-robust standard errors by dyad#
vcov1a <- vcovCL(mod1a, cluster=panel_pres$dyad)
vcov2a <- vcovCL(mod2a, cluster=panel_pres$dyad)
vcov3a <- vcovCL(mod3a, cluster=panel_pres$dyad)
vcov4a <- vcovCL(mod4a, cluster=panel_pres$dyad)

ses <- list(vcov1a, vcov2a, vcov3a, vcov4a)
ses <- lapply(ses, diag)
ses <- lapply(ses, sqrt)

stargazer(mod1a, mod2a, mod3a, mod4a, se=ses,
          keep=c("ideo_prox", "avg_pres_dist", "emp_diff", "avg_house_dist", "total_rules_sum",
                 "overlapping_laws", "total_mentions", "avg_court_dist"),
          covariate.labels = c("Agency Alignment", "Presidential Misalignment",
                               "Court Misalignment", "Overlapping Laws",
                               "Presidential Attention", "Log(Total Rules)",
                               "Employment Difference",
                               "House Misalignment", "Agency Alignment $\\times$ Pres. Misalignment"),
          dep.var.labels.include = T, keep.stat=c("adj.rsq", "n"), digits=3,
          star.cutoffs = c(.05, .01, .001), type = "text")

#Table C3#

mod1y <- lm(weight_binary ~ ideo_prox + avg_pres_dist + dyad + factor(year), data=panel_net)
mod2y <- lm(weight_binary ~ ideo_prox*avg_pres_dist + + dyad + factor(year), data=panel_net)
mod3y <- lm(weight_binary  ~ ideo_prox + avg_pres_dist + avg_court_dist + overlapping_laws + log(total_mentions+1) + log(total_rules_sum) + emp_diff + polit_diff + total_polit + avg_house_dist + dyad + factor(year), data=panel_net)
mod4y <- lm(weight_binary  ~ ideo_prox*avg_pres_dist + avg_court_dist + overlapping_laws +  log(total_mentions+1) + log(total_rules_sum) + emp_diff + polit_diff + total_polit + avg_house_dist + dyad + factor(year), data=panel_net)


#Caluculte cluster-robust standard errors by dyad#
vcov1y <- vcovCL(mod1y, cluster=panel_net$dyad)
vcov2y <- vcovCL(mod2y, cluster=panel_net$dyad)
vcov3y <- vcovCL(mod3y, cluster=panel_net$dyad)
vcov4y <- vcovCL(mod4y, cluster=panel_net$dyad)

ses <- list(vcov1y, vcov2y, vcov3y, vcov4y)
ses <- lapply(ses, diag)
ses <- lapply(ses, sqrt)

stargazer(mod1y, mod2y, mod3y, mod4y, se=ses,
          keep=c("ideo_prox", "avg_pres_dist", "emp_diff", "polit_diff", "avg_house_dist", "total_rules_sum",
                 "overlapping_laws", "total_mentions", "avg_court_dist", "total_polit"),
          dep.var.labels.include = T, keep.stat=c("adj.rsq", "n"), digits=3,
          star.cutoffs = c(.05, .01, .001), type = "text")

#Table C4#
mod1c <- lm(I(log(count_rules+1)) ~ ideo_prox + avg_pres_dist + dyad + pres_term, data=panel_pres_slim)
mod2c <- lm(I(log(count_rules+1))  ~ ideo_prox*avg_pres_dist + + dyad + pres_term, data=panel_pres_slim)
mod3c <- lm(I(log(count_rules+1))  ~ ideo_prox + avg_pres_dist + avg_court_dist + overlapping_laws + log(total_mentions+1) + log(total_rules_sum) + emp_diff + polit_diff + total_polit + avg_house_dist + dyad + pres_term, data=panel_pres_slim)
mod4c <- lm(I(log(count_rules+1))  ~ ideo_prox*avg_pres_dist + avg_court_dist + overlapping_laws + log(total_mentions+1) + log(total_rules_sum) +emp_diff + polit_diff + total_polit + avg_house_dist + dyad + pres_term, data=panel_pres_slim)


#Caluculte cluster-robust standard errors by dyad#
vcov1c <- vcovCL(mod1c, cluster=panel_pres_slim$dyad)
vcov2c <- vcovCL(mod2c, cluster=panel_pres_slim$dyad)
vcov3c <- vcovCL(mod3c, cluster=panel_pres_slim$dyad)
vcov4c <- vcovCL(mod4c, cluster=panel_pres_slim$dyad)

ses <- list(vcov1c, vcov2c, vcov3c, vcov4c)
ses <- lapply(ses, diag)
ses <- lapply(ses, sqrt)

stargazer(mod1c, mod2c, mod3c, mod4c, se=ses,
          keep=c("ideo_prox", "avg_pres_dist", "emp_diff", "polit_diff", "avg_house_dist", "total_rules_sum",
                 "overlapping_laws", "total_mentions", "avg_court_dist", "total_polit"),
          dep.var.labels.include = T, keep.stat=c("adj.rsq", "n"), digits=3,
          star.cutoffs = c(.05, .01, .001), type = "text")

#Table C5#
mod1s <- lm(rule ~ ideo_prox + avg_pres_dist + dyad + pres_term, data=panel_pres_slim[panel_pres_slim$overlapping_laws>0,])
mod2s <- lm(rule ~ ideo_prox*avg_pres_dist + dyad + pres_term, data=panel_pres_slim[panel_pres_slim$overlapping_laws>0,])
mod3s <- lm(rule ~ ideo_prox + avg_pres_dist + avg_court_dist + overlapping_laws + log(total_mentions+1) + log(total_rules_sum) + emp_diff + polit_diff + total_polit + avg_house_dist + dyad + pres_term, data=panel_pres_slim[panel_pres_slim$overlapping_laws>0,])
mod4s <- lm(rule ~ ideo_prox*avg_pres_dist + avg_court_dist + overlapping_laws  + log(total_mentions+1) + log(total_rules_sum) + emp_diff + polit_diff + total_polit + avg_house_dist + dyad + pres_term, data=panel_pres_slim[panel_pres_slim$overlapping_laws>0,])


#Caluculte cluster-robust standard errors by dyad#
vcov1s <- vcovCL(mod1s, cluster=panel_pres_slim$dyad[panel_pres_slim$overlapping_laws>0])
vcov2s <- vcovCL(mod2s, cluster=panel_pres_slim$dyad[panel_pres_slim$overlapping_laws>0])
vcov3s <- vcovCL(mod3s, cluster=panel_pres_slim$dyad[panel_pres_slim$overlapping_laws>0])
vcov4s <- vcovCL(mod4s, cluster=panel_pres_slim$dyad[panel_pres_slim$overlapping_laws>0])

ses <- list(vcov1s, vcov2s, vcov3s, vcov4s)
ses <- lapply(ses, diag)
ses <- lapply(ses, sqrt)

stargazer(mod1s, mod2s, mod3s, mod4s, se=ses,
          keep=c("ideo_prox", "avg_pres_dist", "emp_diff", "polit_diff", "avg_house_dist", "total_rules_sum",
                 "overlapping_laws", "total_mentions", "avg_court_dist", "total_polit"),
          star.cutoffs = c(.05, .01, .001), type = "text")

#Table C6#
mod1s <- lm(rule ~ ideo_prox + min_pres_dist + dyad + pres_term, data=panel_pres_slim)
mod2s <- lm(rule ~ ideo_prox*min_pres_dist + dyad + pres_term, data=panel_pres_slim)
mod3s <- lm(rule ~ ideo_prox + min_pres_dist + avg_court_dist + overlapping_laws + log(total_mentions+1) + log(total_rules_sum) + emp_diff + polit_diff + avg_house_dist + dyad + pres_term, data=panel_pres_slim)
mod4s <- lm(rule ~ ideo_prox*min_pres_dist + avg_court_dist + overlapping_laws  + log(total_mentions+1) + log(total_rules_sum) + emp_diff + polit_diff + avg_house_dist + dyad + pres_term, data=panel_pres_slim)


#Caluculte cluster-robust standard errors by dyad#
vcov1s <- vcovCL(mod1s, cluster=panel_pres_slim$dyad)
vcov2s <- vcovCL(mod2s, cluster=panel_pres_slim$dyad)
vcov3s <- vcovCL(mod3s, cluster=panel_pres_slim$dyad)
vcov4s <- vcovCL(mod4s, cluster=panel_pres_slim$dyad)

ses <- list(vcov1s, vcov2s, vcov3s, vcov4s)
ses <- lapply(ses, diag)
ses <- lapply(ses, sqrt)

stargazer(mod1s, mod2s, mod3s, mod4s, se=ses,
          keep=c("ideo_prox", "min_pres_dist", "emp_diff", "polit_diff", "avg_house_dist", "total_rules_sum",
                 "overlapping_laws", "total_mentions", "avg_court_dist"),
          covariate.labels = c("Agency Alignment", "Presidential Misalignment",
                               "Court Misalignment", "Overlapping Laws",
                               "Presidential Attention", "Log(Total Rules)",
                               "Employment Difference", "Politicization Difference",
                               "House Misalignment", "Agency Alignment $\\times$ Pres. Misalignment"),
          dep.var.labels.include = T, keep.stat=c("adj.rsq", "n"), digits=3,
          star.cutoffs = c(.05, .01, .001), type = "text")

#Figure C1#

library(lfe)

mods <- as.list(rep(NA, length(unique(panel_pres_slim$dyad))))

#Estimate main model dropping one dyad at a time#
for(i in 1:length(unique(panel_pres_slim$dyad))){
  mods[[i]] <- felm(rule ~ ideo_prox*avg_pres_dist + avg_court_dist + overlapping_laws + 
                      log(total_mentions+1) + log(total_rules_sum) + emp_diff +
                      polit_diff + total_polit + avg_house_dist | dyad + pres_term |0| dyad,
                    data=subset(panel_pres_slim, dyad!=unique(panel_pres_slim$dyad)[i]))
  message(i,"/",length(unique(panel_pres_slim$dyad)))
}

coefs <- lapply(mods, coef)

int <- NA
for(i in 1:length(coefs)){
  int[i] <- coefs[[i]][11]
}

ints <- data.frame(int=int)
ints$se <- NA
for(i in 1:length(coefs)){
  ints$se[i] <- mods[[i]]$se[11]
}


ggplot(ints) + geom_point(aes(x=1:496, y=int)) + 
  geom_errorbar(aes(x=1:496, ymin=int-1.96*se, ymax=int+1.96*se)) +
  ylim(c(0,2.5)) + theme_bw() + xlab("") + theme(axis.text.x = element_blank()) +
  ylab("Estimated Coefficient") + geom_hline(yintercept=0, lty=2, size=1, color="red") +
  ggtitle("Dropping One Dyad at a Time")


#Plot distribution of coefficients#

ggplot(ints) + geom_density(aes(x=int), fill="grey", alpha=.5) + theme_bw() +
  geom_vline(xintercept=coef(mod4)[length(coef(mod4))], lwd=1) + xlab("Estimated Coefficient") +
  ylab(" ") + xlim(c(0, 2)) + ggtitle("Dropping One Dyad at a Time")  +
  theme(panel.grid = element_blank()) + geom_vline(xintercept=0, color="red", lty=2) +
  geom_vline(xintercept = mean(ints$int), lty=2)


mods2 <- as.list(rep(NA, 32))
agencies <- unique(c(panel_pres_slim$agency_1,
                     panel_pres_slim$agency_2))
for(i in 1:length(agencies)){
  mods2[[i]] <- felm(rule ~ ideo_prox*avg_pres_dist + avg_court_dist + overlapping_laws +
                       log(total_mentions+1) +log(total_rules_sum) + total_polit + emp_diff +
                       polit_diff + avg_house_dist | dyad + pres_term |0| dyad,
                     data=subset(panel_pres_slim, agency_1!=agencies[i] & agency_2!=agencies[i]))
  message(i, "/", length(unique(c(panel_pres_slim$agency_1, panel_pres_slim$agency_2))))
}

coefs2 <- lapply(mods2, coef)

int2 <- NA
for(i in 1:32){
  int2[i] <- coefs2[[i]][11]
}

ints2 <- data.frame(int=int2)

ints2$se <- NA
for(i in 1:length(coefs2)){
  ints2$se[i] <- mods2[[i]]$se[11]
}

ggplot(ints2) + geom_point(aes(x=1:32, y=int)) + 
  geom_errorbar(aes(x=1:32, ymin=int-1.96*se, ymax=int+1.96*se)) +
  ylim(c(0,2.5)) + theme_bw() + xlab("") + theme(axis.text.x = element_blank()) +
  ylab("Estimated Coefficient") + geom_hline(yintercept=0, lty=2, size=1, color="red") +
  ggtitle("Dropping One Agency at a Time")


ggplot(ints2) + geom_density(aes(x=int), fill="grey", alpha=.5) + theme_bw() +
  geom_vline(xintercept=coef(mod4)[length(coef(mod4))], lwd=1) + xlab("Estimated Coefficient") +
  ylab(" ") + xlim(c(0, 2)) + ggtitle("Dropping One Agency at a Time")  +
  theme(panel.grid = element_blank())  + geom_vline(xintercept=0, color="red", lty=2) +
  geom_vline(xintercept = mean(ints2$int), lty=2)

#Figure D1#
FEadj <- data.frame(std=scale(panel_pres_slim$ideo_prox_raw, center=T, scale=F))
FEadj$FEadj <- lm(ideo_prox_raw ~ dyad + pres_term, data=panel_pres_slim)$residuals

ggplot(FEadj) + geom_density(aes(x=std), fill="grey", alpha=.5) +
  geom_density(aes(x=FEadj), fill="red", alpha=.5) + theme_bw() +
  ggtitle("Agency Alignment (Centered on Zero)") + ylab("Density") +
  annotate(-.25, 3.75, geom="text", label="After FEs\nSD=0.103", col="red") + xlab(" ") +
  annotate(.5, 3.75, geom="text", label="Original Distribution\nSD=0.135") + xlim(c(-1,1)) +
  theme(panel.grid = element_blank())


FEadj_pres <- data.frame(std=scale(panel_pres_slim$avg_pres_dist_raw, center=T, scale=F))
FEadj_pres$FEadj <- lm(avg_pres_dist_raw ~ dyad + pres_term, data=panel_pres_slim)$residuals

ggplot(FEadj_pres) + geom_density(aes(x=std), fill="grey", alpha=.5) +
  geom_density(aes(x=FEadj), fill="red", alpha=.5) + theme_bw() +
  ggtitle("Presidential Misalignment (centered on Zero)") + ylab("Density") +
  annotate(.15,4,geom="text", label="After FEs\nSD=0.087", col="red") + xlab(" ") +
  annotate(-.35, 2, geom="text", label="Original Distribution\nSD=0.307") + xlim(c(-.6,.6)) +
  theme(panel.grid = element_blank())


#Figure D2#
within <- data.frame(within=tapply(panel_pres_slim$ideo_prox_raw, panel_pres_slim$dyad, max)-tapply(panel_pres_slim$ideo_prox_raw, panel_pres_slim$dyad, min))
ggplot(within) + geom_histogram(aes(x=within), fill="grey", bins=20, col="black") +
  theme_bw() + geom_vline(xintercept = median(within$within), lty=2) + ylim(c(0,200)) +
  geom_vline(xintercept = quantile(within$within, .95), lty=2) +
  annotate(.3, 186, geom="text", label="Median") + ylab("Frequency") +
  ggtitle("Within-Dyad Range of Agency Alignment") + xlab(" ") +
  annotate(.6, 105, geom="text", label="95th\nPercentile") +
  geom_segment(aes(x=.3, xend=median(within)+.01, y=179,yend=179),
               arrow=arrow(angle=30, ends = "last", length=unit(.1, "in"))) +
  geom_segment(aes(x=.6, xend=quantile(within,.95)+.01, y=90,yend=90),
               arrow=arrow(angle=30, ends = "last", length=unit(.1, "in"))) +
  theme(panel.grid = element_blank())


within_pres <- data.frame(within=tapply(panel_pres_slim$avg_pres_dist_raw, panel_pres_slim$dyad, max)-tapply(panel_pres_slim$avg_pres_dist_raw, panel_pres_slim$dyad, min))
ggplot(within_pres) + geom_histogram(aes(x=within), fill="grey", bins=20, col="black") +
  theme_bw() + geom_vline(xintercept=median(within_pres$within), lty=2)+ ylim(c(0,250))+
  annotate(.8, 240, geom="text", label="Median") + ylab("Count") +
  ggtitle("Within-Dyad Range of Presidential Misalignment") + xlab(" ") +
  geom_vline(xintercept=quantile(within_pres$within, .95), lty=2) +
  annotate(.97, 130, geom="text", label="95th\nPercentile") +
  geom_segment(aes(x=.83, xend=median(within)+.01, 233, yend=233),
               arrow=arrow(angle=30, ends="last", length=unit(.1, "in"))) +
  geom_segment(aes(x=1, xend=quantile(within, .95)+.01, y=112, yend=112),
               arrow=arrow(angle=30, ends="last", length=unit(.1, "in"))) +
  theme(panel.grid = element_blank())

#Figure E1#
panel_net <- subset(panel, dyad %in% panel_pres_slim$dyad)
net <- data.frame(weight = tapply(panel_net$weight, as.character(panel_net$dyad), sum))

#Join dyad-level variables#
net$dyad <- rownames(net)
net <- join(net, panel_net, by="dyad", type="left", match="first")
net$agency_1 <- as.character(net$agency_1)
net$agency_2 <- as.character(net$agency_2)
net2 <- net[,c(4,5,2)]
net2 <- subset(net2, weight>0)

#Create undirected network object#
g <- graph_from_data_frame(net2, directed=F)

#Color agency types#
V(g)$color <- grepl("department", vertex_attr(g)[[1]], ignore.case = T) 
V(g)$color <- ifelse(V(g)$color==TRUE, "black", "white")


#Plot network -- Graphing has a random component that cannot be controlled with seeds#
#so this is commented out to not overwrite the clearest plot#

plot(g, layout=layout_with_lgl, vertex.label=NA,
      vertex.size=degree(g)/3, edge.width=log(E(g)$weight))
 title(main="Pooled", cex.main=2, adj=0)

#Create yearly edge lists#
netsTerm <- as.list(rep(NA, 4))
netsTerm2 <- as.list(rep(NA, 4))

for(i in 1:length(unique(panel_pres_slim$pres_term))){
  netsTerm[[i]] <- subset(panel_pres_slim, pres_term==unique(panel_pres_slim$pres_term)[i])
  netsTerm[[i]] <- netsTerm[[i]][,c(2,3,19)]
  netsTerm2[[i]] <- subset(netsTerm[[i]], count_rules>0)
  
}

#Create network objects for each year#
gBushI <- graph_from_data_frame(netsTerm2[[1]], directed=F)
gBushII <- graph_from_data_frame(netsTerm2[[2]], directed=F)
gClintonII <- graph_from_data_frame(netsTerm2[[3]], directed=F)
gObamaI <- graph_from_data_frame(netsTerm2[[4]], directed=F)

#Color nodes#
for(i in 1:length(unique(panel_pres_slim$pres_term))){
  gTerm <- get(gsub("[[:space:]]", "",paste0("g", unique(panel_pres_slim$pres_term)[i], sep="")))
  V(gTerm)$color <- grepl("department", vertex_attr(gTerm)[[1]], ignore.case = T) #Color cabinet departments black
  V(gTerm)$color <- ifelse(V(gTerm)$color==TRUE, "black", "white")
  V(gTerm)$color2 <- grepl("office", vertex_attr(gTerm)[[1]], ignore.case = T)
  V(gTerm)$color <- ifelse(V(gTerm)$color2==TRUE, "grey", V(gTerm)$color)
  assign(gsub("[[:space:]]", "",paste0("g", unique(panel_pres_slim$pres_term)[i], sep="")), gTerm)
}

#Plot network over time #
par(mar=c(1.5,0.5,2.5,0.5),mfrow=c(1,4))
for(i in 1:length(c("Clinton II", "Bush I", "Bush II", "Obama I"))){
  plot(get(gsub("[[:space:]]", "",paste0("g", c("ClintonII", "BushI", "BushII", "ObamaI")[i], sep=""))), vertex.label=NA,
       vertex.size=degree(get(gsub("[[:space:]]", "",paste0("g", c("ClintonII", "BushI", "BushII", "ObamaI")[i], sep=""))))/1.5,
       edge.width=log(E(get(gsub("[[:space:]]", "",paste0("g", c("ClintonII", "BushI", "BushII", "ObamaI")[i], sep=""))))$count_rules), frame=T)
  title(c("Clinton II", "Bush I", "Bush II", "Obama I")[i], adj=0, line=.5, cex.main=2.5)
}
dev.off()

